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ABSTRACT 


This thesis is concerned with the analysis of various methods for the numeri- 
cal solution of the shallow water equations along with the stability of these methods. 
Most of the thesis is concerned with the background and formulation of the shallow 
water equations. The derivation of the basic equations will be given, in the primitive 
variable and vorticity-divergence formulation. Also the shallow water equations will 
be written in spherical coordinates. Two main types of methods used in approximat- 
ing differential equations of this nature will be discussed. The two schemes are finite 
difference method (FDM) and the finite element method (FEM). After presenting 
the shallow water equations in several formulations, some examples will be presented. 
The use of the Fourier transform to find the solution of a semidiscrete analog of the 


shallow water equations is also demonstrated. 
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I. INTRODUCTION 


The majority of this work is a compilation of past mathematical papers con- 
cerning the numerical solution of the shallow water equations, and I can in no way 
take credit for them. Hopefully I have put this material in a readable and under- 
standable text that will provide a good basis for others to use for reference in this 
area of mathematics. The next several paragraphs will highlight the key points in 
slightly more detail than the Table of Contents. 

Chapter II shall review and discuss the mathematical formulation of flow 
in shallow regions. Two formulations will be given, namely the primitive and the 
vorticity-divergence formulations of the shallow water equations. It will also look at 
the formulation of the equations in spherical coordinates. 

Chapter III will look at the linearization of the shallow water equations. In 
Cartesian coordinates we discuss the one dimensional and two dimensional cases. 
The linearized vorticity-divergence form of the shallow water equations is derived in 
this chapter for the two dimensional case. The spherical case is also discussed here. 
Note that in the Cartesian case, the linearized equations have constant coefficients in 
contrast to the spherical coordinates. This means that linear stability analysis will 
be more difficult in the latter. 

Chapter IV shall develop two of the different types of approximations com- 
monly used. The first is the finite difference, and the second is the finite element. The 
finite element may also be referenced to as the Galerkin method and will be referred 
to by both names in this thesis. Other families of methods such as spectral or finite 
volume approximation are beyond the scope of this thesis. 

In Chapter V, the stability analysis of the shallow water equations will be 
examined. For this analysis we use the linearized equations obtained in Chapter II. 
Fourier transforms in a one dimensional case will be covered. Fourier transforms 


will also be used to look at the two dimensional case. The spherical case requires 





special consideration as in Longuet-Higgins(see [Ref. 1]) or Neta (see [Ref. 2]). A 
summary of the results for the various techniques are provided in a table for ease of 
later examination. 

Also covered in this report will be an extensive bibliography and reference 
section, so that future students studying in this area will be able to pick up more 


easily where this thesis leaves off. 








Il. SHALLOW WATER MODEL 


A. MODEL BACKGROUND 

For this model consider a sheet of fluid with constant and uniform density. 
(See [Ref. 3].) The height of the surface of the fluid above the reference level z = 0 is 
h(x, y,t). We model the body force arising from the potential ¢ = gh with atmosphere 
or ocean in mind. ¢ is a vector directed perpendicular to the z = 0 surface, or db can 
be said to be antiparallel to the vertical axis (i.e., db is in the direction opposite to 
the vertical axis). The rotation axis of the fluid is the z-axis in this model. In this 
case the Coriolis parameter f is 20 since @ = kN. The rigid bottom 1s defined by 
the surface z = hg(x,y). The velocity has components u, v, and w in the z, y, and 
z directions respectively. Though the pressure of the fluid surface can be arbitrarily 
imposed, for this model it will be assumed to be constant. Lastly, the fluid is assumed 
inviscid, in other words, only the motions for which viscosity is not important are 
considered. 

In this model, because the depth of the fluid, h — hg, varies over time or space, 
let H be the average depth of the fluid. H characterizes the vertical scale of motion 
also. Let L be the characteristic horizontal scale for the motion. Then a fundamental 
condition which will characterize shallow water theory will be 


= <1 
L 


which is also called the hydrostatic approximation with long wavelengths. (See [Ref. 
4].) The shallow water model contains several of the important dynamical features 
of the atmosphere and ocean while being simple enough to be easily understood. 
The major physical difference with this model and reality is the absence of density 
stratification that is present in the real fluids such as Earth’s atmosphere or oceans. 
The hydrostatic approximation also allows p to vary with z, but we will consider p 


constant for this model. 














Recalling the equation for motions of fluids from Haltiner and Williams (see 
[Ref. 5], we have 


dV a ; 
<= -a¥p— 2x V4 g+h (2.1) 


a is the specific gravity (i.e. :)> g is the sum of gravitational and centrifugal forces, 
and F is the force due to friction. For our model we will be neglecting g and F 
because with our assumptions they are much smaller in magnitude than the Coriolis 


forces. Recall that 
V(V-V)+(VxV)xV, (2.2) 


where V = ua + vj + wk. Using these two basic equations, we will construct our 


model, 


B. MODEL EQUATIONS 


There are several formulations in the literature. We will cover 


® primitive variable formulation 


@ vorticity divergence formulation. 
The formulations using stream function and velocity potential will not be discussed 
here, (see, for example, [Ref. 5]). 
1. Primitive Form 


Now follow the consequences of the model in the realm of the dynamical equa- 
tions of motion. Recalling the Navier-Stokes equations which describe the conser- 


vation of mass and momentum (see [Ref. 6]), we see that the dynamics and the 


thermodynamics decouple due to the specification of incompressibility and constant 

















density. (See [Ref. 3].) This reduces the equation of mass conservation to the condi- 
tion of incompressibility: 

ot get Geo. (2.3) 
Now the first two terms of this equation are of order . where U can be considered 
the characteristic scale for the horizontal velocity. It then follows that the scale for 


the vertical velocity (W) is smaller than or equal to the order dU. This represents an 


upper bound for the vertical velocity, and it can be smaller than order 6U if there is 


Ou Ov 
cancellation between — and —. 
Ox Oy 
Now an estimate of the momentum equations in component form are: 
Ou Ou Ou Ou 1 Op 


— +u— +0—+u——-—fv = 


At ' Oz Oy Oz pp Ox 

yu ue ou uw 4 P 

rT A L 6A pL 

OU eye OY fa Oe — _1 9p 

a Oe Oy Oz op Oy a 
Uw UW 4 P 7 
ar Fr i w oL 

Ow | Ow aa _ _1 Op 

Ot Ox Oy Oz pp Oz 

w Uw UW WwW" a2 

T L L H pH 


f will be defined differently depending on the form of the shallow water equations. 
In (2.4) each term has the order of magnitude written immediately below it in terms 
of the characteristic scales where J’ is the characteristic scale for time and P is the 
characteristic scale for the pressure field. Equation (2.4) follows from (2.1) and (2.2) 


in the following manner. Remember that we are neglecting the gravity term and the 





force due to friction. Now when we expand out (2.1), we get 


dV 1 Op- 1 Op- 1 Op; 
di p Ox poy” p Oz 


which equals 


(2.5) 
+(2w cos y — fv)it fuj — 2uQ-cos ok. 

Since we are modeling a rectangular system first, the y terms drop out and (2.5) 
becomes . 

dV 1 Op-+ 1 Op+ 1 Op hd cg zs 

a ey — eK : 2.6 

dt p Ox p dy” poz Just fuy (2.6) 
When we look at the spherical form of the shallow water equations, the y terms are 


used. Note that 


Now, expanding out (2.2) we get 


dV _ OV (Ou, du, Ow\s, (du, dv, dw) - 
"ax Ox Ox . ay 


a OL 
(ed 4 v2 4 wd!) gy fy (Ou Ow) _, (8 _ Ou] = 
"Oz dz Oz "\Oz 7 Oz \ as - Oy 


in «(3S -F Hi ee lee a ee 
dr dy dy dz}i7 "| \ay az z Oxf] 


Collecting terms, this becomes 





or more simply 


dV OV = 
eo ee 


Equating the right hand side of (2.6) and this equation and then separating into the 


three components will give us (2.4). 


Note that the total pressure (p) is 


p(x, y,2,t) = —pgz + p(x, y, z,t). 


The horizontal pressure gradient is independent of z if one uses the hydrostatic ap- 
proximation 
Op 


Pian as O(5") (2.7) 


H es ; 
where 9¢ is 7 and 6 <1. This approximation follows from the scale analysis of 
the momentum equations and the incompressibility condition. Integrating (2.7) and 


using the boundary condition 
p(t, y,h) = po 
_ yields 
p = pg(h — z) + po. 


A way to look at this is to think of the pressure which is in excess of pp at any point 
simply as the weight of the unit column of fluid above that point at that particular 
instant in time. Remember that the horizontal pressure gradient is independent of z, 


and therefore the horizontal accelerations are independent of z. 


U 


The Rossby-wave number is the ratio of inertia to the Coriolis terms or (ap) 
Recalling the Taylor-Proudman theorem (see, [Ref. 3] p. 43) which simply put states 


that if the Rossby-wave number is small, friction can be ignored, and baroclinic vector 
is identically zero (i.e., there is no additional pressure change with a change in height 


as in this model), then it follows that the horizontal accelerations are independent of 


"7 
Cd 


Ou Ov 


,, ama ae 





In this model, the Taylor-Proudman theorem applied to a homogeneous fluid requires 


the velocities to be independent of z. This allows the horizontal momentum equations 


to become 


Ow Ou Oe — _1 Op 
dt ace ay °° ~~ pdx 
OU goa. OF 5 —_ _1 Op 
bt Or Oy ~ 5 Ay 


(2.8) 


Let us look at the scale analysis of the vertical momentum equation of (2.4). 


Consider w to be of O (#) = O(6). Scale u and v by U, and w by 


know that 
ai | | Oi 
Ox! Oy’ Oz’! 


where 


ss Zz 
w = 6Uw Ze = —. 
a 
The vertical momentum equation can be written as 


pf os ne ee ap Ow 
Ot L Oz! L Oy' H Oz! 


Multiply through by H and this becomes 





Ow Ow Ow =. _ 


6U WS + Uoun + Se, + wr 


#U. By (2.3) we 


11 Op 
«Ap Oz" 


_1 op 
pp Ozh 


Removing all the terms of O(5?), note from (2.7) that p is of O(67) also, we have 


Ow 
bUH ~~ = 0. 


Therefore the vertical momentum equation is an identity. 


Now we will rewrite the incompressibilty condition (2.3), 


bw (du, dv 
Oz Or Oy) 


Since the horizontal accelerations are independent of z, we can integrate the above 
relationship to get 


Ou Ov : 
w(z,y, Z,t) se t + = + w(z,y,t). (2.9) 


The condition of no normal flow at the bottom requires 


It naturally follows that 


Ou dv Ohg Ohp 
Oy 


(wWx,y,2,t) = (hp — z) (sa a "Ay 


dz 
The vertical velocity w(x, y, z,t) = — at the upper surface z = h represents the rate 


dt 
at which the free surface is rising. Thus w(z,y,h,t) = oe and (2.9) become 
Oh Oh Oh 
t)= — aera a 
w(zy,h,t) a “Oe Oy 
Therefore 
Oh O 0 
— [u(h—h — ~~ = 
ae + a [u(h— ha)] + [o(h ~ he} =0 


This equation, combined with the horizontal momentum equations (2.9) form the 


shallow water equations. 


Ov Ov Ov dh 

Ot Or "ay ne = og By (2.10) 
Oh O a 

OL Bs Or [u(h — hg)| = ay [o(h _ hp)| = (0. 


u, v, and fh are called the primitive variables and (2.10) is the primitive form of 


the shallow water equations. 








2.  Vorticity-divergence Form 


For vorticity-divergence, we start with the definitions. The relative vorticity, 


denoted by ¢, is defined by 
C = Ug — Uy, (2.11) 


and the absolute vorticity, denoted by Q, is given by 
Q=c+f. (2.12) 
The divergence, denoted by D, is given by 
DP Ve. (2.13) 


and the kinetic energy (per unit mass), denoted by K, is given by 


u? + v? 
2? 





K= (2.14) 


To get the vorticity-divergence formulation, we start by differentiating the first 
equation of (2.10) with respect to y and subtract it from the second equation of (2.10) 


differentiated with respect to x gives 


Vet bt Uses + UVez + Ugvy al UUzy a fue i Ghy (2 15) 


—Uyt — Uy, — UUgy — Vylly — VUyy + foy — ghey = 0. 
Rearranging and cancelling terms will yield 


(vz — Uy), + Uz(Vz — Uy) + u(vz — ty), 


+v,(vz — Uy) + v(vz — Uy), + fluc + vy) = 0. 


Rewriting in terms of (2.11) and (2.13) gives 
G+urCtul,+v¢C+vG+ fD=0. 
This can be rewritten in the form 
Gr + (uC), + (vd), + fD = 0 


10 





or 
ac 
— D=0 
at! 
This is called the vorticity equation. To get the divergence equation, we difffer- 
entiate the first equation of (2.10) with respect to z and add to it the second equation 


of (2.10) differentiated with respect to y gives 


Ust + ist + Wee + Ugly + UUzry — fvz = hex 


(2.16) 
tye + UyVe + UVay + Vy? + UVyy + fully + ghyy = 0. 


Again, rearranging terms will yield 


(Uz tr Vy), a U( Us i Vy). ag v(Uz a Vy) + 
Ug? + Vy? + 2uyVe — f(ve — Uy) + (hee + hy) = 0. 


Using (2.11) and (2.13) this can be further reduced to 
D,+ uD, + vDy4+ Viet Vy" + Quyvz -— fC + gV*h =*(); 


The vorticity-divergence form of the shallow water equations in terms of rela- 


tive vorticity and divergence are 


Cr + (uC), + (v¢), + FD . = 0 
D,+uDz + uD, + uz? + vy? + Quyve-fC+gV?h = 0 (2.17) 
ht + [u(h — he), + [v(h — he), = 


To simplify a bit further, we can use (2.11) through (2.14) on (2.15) and (2.13). 


First, rearrange and cancel terms of (2.15) in the following manner. 


Urt — Uyt oF UyUr — UyUy + iu; i UVgr — Utyrt 


VyVz — Vyty + fly + VV, — VUy, = 0 


Rearranging further gives 


(vz — uy), + Uc(ve — Uy + f) + Ulver — Uyr)+ 
Vy(Vz — Uy + f) + v(vzy — Uyy) = 0. 


1] 





This leads to 
G +ucQ + uQz t+ vyQ + vQ, = 0, 


OT 


Cr + (uQ), + (vQ), = 0. (2.18) 


This can also be written as 
Q: + (uQ), + (v@), = 0 


since f; = 0. 
Now we simplify (2.16) in a similar fashion. First, rewrite the equation in the 


following manner. 


Unt + Vyt + Ghee + Ghyy + Utes + Ugle + UUVz_ + VeVet 


Notice that we added in and subtracted some terms that were alike. Now reduce this 


equation to 
Urt + Vyt + Ghee + ghyy + (wus + VUz)_+ 
(uty + Vvy), — Use + Vgty — fe — VU eet 
VUye + UyVe — Uytly + fly + UVgy — UUyy = 0. 


Reorganize the terms to further reduce to 


u? + vy? 
Unt = Vyt ai Ghixs = Ghyy =e (=F) + 


Z 
u? + y? 
( ; — v,(vz — uy + f) — v(vz — uy + f),.+ 
~ yy 





uy(vr — Uy + f) + ulvz —uy+ f), = 0. 


Now use the definitions of Q and Av to reduce this equation to 


(ur + vy), + g(her + hy) + Nee + Ky 
vzQ — vQ, + uyQ + uQ, = 0. 


12 











Thus our equation becomes 


Dit+gV7h+V7K — (vQ)2 + (uQ), = 0. (2.19) 


Combining (2.18), (2.19), and the last equation of (2.10) we have the shallow 


water equations in the vorticity divergence form. 


C+ (uQ), + (vQ), = 10 
Di+gV7h+ V?K —(vQ)2+(uQ)y = 0 (2.20) 
hy + [u(h — hp)], + [v(h — he)], = 0 


C. SPHERICAL COORDINATES 
1. Background 


Where Cartesian coordinates locate a point using 3 vectors, spherical coor- 
dinates locate the same point using two angles and a distance. Typically, spherical 
coordinates represent a point in space with an ordered triplet such as (a, ¢, A). These 


variables are related to (x,y,z) by 


x= asin ¢cos A 
y = asin dsin A — (2.21) 
z= acos@. 
Let us consider a channel that goes around the whole Earth. In this model, 
let A be the longitude, let ¢ be the latitude, let r be the radial distance, let a be the 


average radius of the Earth, and let z be the average sea level. Now, recall equations 


(2.5) and (2.2). 


+(2w cos ¢ — foyit fuj — 2uf cos ok. 


13 


and _ _ 

dV OV las o a ee 
| ae ap vg a ee 
Using results drawn from Haltiner and Williams (see [Ref. 5]), and neglecting the 


force due to friction, F’, gives us 


fe age butt 
~ ae ae ay" Gz)" 


+ hee oa eae } + Rl hc k 
az}? ot Ox Oy Oz 


Ot Ox Oy 
a POP 29 + —— (vsing —w os d) 1 (2.22) 
p Or acos : ) 





+ {222 + (204 : ) usin o— 2 5 
poy a cos a, 


2 
+ ae 20 + : ucosé + — k; 
p Oz acos @ a 


Separating along the three axes and we have 














Ou Ou Ou Ou 1 0p u 

OS a ae an ee Oe (20 + +) (usin d — wcos ¢) 
Ov Ov Ov Ov 10p u vw 

= use oe gt + ASE (20+ 2) using — ee 
Ou Ow dw Ow  10p u v? 

0m Oe gu go we 4 OF + 0+ (204 A) woos dt 


oa 10 lo : Oh 
Note that p is independent of z, and =~ and - can be written in terms of ox 
pox p Oy Ox 


ol 
and I respectively. Also note that 
Y 


du _ dud , duds 
Or = OAx_-~—ssOH Oz 


Ou ] rece ] 
~ O@—asingsinA Odacosdcos dr 
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and similarly for the other partials. Substituting these and the Coriolis parameter, 
f = 2Msin 8, into (2.23) gives us the shallow water equations in spherical coordinates. 


(See [Ref. 7].) 








Oo ee got —|f+et 5 pe 

Ot acosé es 06 a sea acos@@\ 

Ov 1 Ov Ov U gOh | 

af an So + vcs 052 + [f+ “tano|u+ 25% =o (2.24) 
Oh 1 [a 8 ; 

Di cose (hn) 4- 9g (he COS 2) = 0. 


Numerical approximations of the shallow water equations in spherical coordinates are 
given for example by Turkel and Zwas (see [Ref. 8]) and Arakawa (see [Ref. 9]) et 
al. Analysis of the schemes were given by Neta and Navon (see [Ref. 10]) et al. 
Navon and de Villiers et al. have applied the Turkel-Zwas scheme to a hemispheric 


barotropic model (see [Ref. 11]). 
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IIT. LINEARIZATION 


In this chapter, we will see how to linearize the shallow water equations in 
Cartesian and spherical coordinates. The Cartesian coordinate case starts with either 


the primitive or the vorticity-divergence formulation. 


A. TWO DIMENSIONS 
1. Basic Formula 


First separate each of the three variables into the average which we can consider 


constant plus the perturbation in the variable. This should look like: 


u=Us+u 
v=V4u' (3.1) 
h=H+h’ 


where the capital letters are used for the average or mean values and the primed 
variables are the perturbations. Now substituting this into (2.7) and omitting all the 


second order terms gives 


dul du’ | dw’ O(H + h’) | 

Eases, | a ace = / = 

a! Bp Oy ae ee Dee 0 

Ov, Ov" Ov! O(H +h’) 

jen | Sean t eee f 2 

a tg, tg, t fl tito 0 

on a a. 9 

aan pee — Renee a ] yp, G r = 

a + ar [u’( H ha)l+ > [o'(h—hp)|)+(U +N oe [fH +h he) 0. 
(3.2) 


For simplicity and without loss of generality we sometimes consider a zero 
mean flow, i.e. Uo = V = 0. This and dropping the primes will give the linearized 


form of the shallow water equations which follows. 
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Ov Oh 

at I¥t 9g, =9 

Oh Ou Ov 

3 th ha) (Fe +52) =0 


If topography is not considered, then the linearized shallow water equations 


can be written 


Ou Oh 
Ot — fut I7- =i 
Ov Oh 
ied ee) 
ot cen 7 Oy 
ah, y (Ou, ar) _y 
Ot Ox Oy 
This can also be written as 
ur— fut ghz = iV 
utfutgh, = 0 (3.3) 


h,+ H(uz+vy) = 0. 
Another source for the development of the linearized shallow water equations is Gill. 


(See [Ref. 4]). 


2.  Vorticity Divergence Formula 

The two dimensional linearized vorticity divergence form of the shallow water 
equations is obtained by using (2.8) and (2.9) on (3.3). It should be noted that 
the linearized vorticity divergence formula can be found linearizing (2.16) or (2.19). 
Taking the partial differential of the first equation of (3.3) with respect to y and 
subtracting it from the partial differential of the second equation of (3.3) with respect 
to x yields 

Vie + fuc + ghey — Uy + fly — ghey = 9. 
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With a little algebra this equation can be rewritten as 
(vz — ty), + f (te + vc) = 0. 
which when using the definitions of ¢ and D reduces to 
G+ fD=0. (3.4) 


This is the first equation of the linearized vorticity divergence form. By taking the 
partial of the first equation of (3.3) with respect to x and adding it to the partial 
differential of the second equation of (3.3) with respect to y we get 


Ure — {2 + ghee + Viy + fuy + ghyy = 0. 
With a little algebraic simplification our equation becomes 
(te + vy), — f(v2 — Uy) + ghee + hyy) = 0. 
Again, using the definitions of ¢ and D, this reduces to 
D,—fC+gV*h =0. (3.5) 


(3.4) and (3.5) along with the last equation of (3.3) form the two dimensional lin- 


earized vorticity divergence form of the shallow water equations. 


G+ fD =) 
D,—fC+gqV?h = 0 (3.6) 
hy+ H(uz + vy) = 0 
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B. LINEARIZED SHALLOW WATER EQUATIONS IN 
ONE DIMENSION 


In the one dimensional case, all dependence on y is eliminated and we end up 


with P oh 
U 
a 1° +99, =° 


Ov 


a + fu=0 
or more simply 
u:—fvtgh, = 0 
v,+ fu = 0 (3.7) 
hit Hu, =), 


C. SPHERICAL 
To linearize the spherical version of the shallow water equations, simply take 
(3.1) and substitute it into (2.26). The resulting equation follows quite easily once 


all second and higher order terms are omitted. 


Ou g Oh 





ot arr. ala 

Ov g Oh 

Ry tiut ia = 

Oh H Ou O 

Or Geos (= + Syl cos ) = (). 


Notice that the linearized equations have nonconstant coefficients. This complicates 


the linear stabilty analysis. See Longuet-Higgins [Ref. 1] and Neta [Ref. 12] et al. 
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IV. APPROXIMATIONS 


A. FINITE DIFFERENCES 
1. Introduction 


One of the first steps in using finite difference methods is to replace the con- 
tinuous problem domain by a difference mesh or a grid. Let f(x) be a function of 
the single independent variable x for a < x < b. The interval [a, }] is discretized by 
considering the nodes a = rg < 21 < +++ < EN < Zn41 = 6, and we denote f(x;) by 
f;. The mesh size is 234, — 2; and we shall assume for simplicity that the mesh size 


is a constant 


h= 





N+1 


and 


e;=a+th i=0,1,---,N+41 


In the two dimensional case, the function f(z,y) may be specified at nodal 
point (z;,y;) by fi;. The spacing in the x direction is h, and in the y direction is hy. 
Taylor series expansions of functions in several variables play an important 
role in formulation and classification of finite difference methods. The Taylor series 


expansion for f;,; about the point z; is given by 
h? he 
Fe Teta Gh Pade Oe 


The Taylor series expansion for f;41;4; about the point (2z;,y;) is given by 


h? h? 
fizijar = Sig + (Arte + hy fy)iz + (> See + hehy fey + 5 Suv) 4... 


2. Finite Differences 
Using text drawn from Neta (see [Ref. 13]), we will cover certain aspects of 
finite differencing. An infinite number of difference representations can be found for - 


the partial derivatives of f(z,y). Let us use the following operators: 


AcSis = Siary — Si 1* forward difference operator 
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Veli; = fij — fi-1;. 1° backward difference operator 
Oefij = fittg — fi-1j centered difference 
Orbis = Ferasa3 — Si-1/2; 
Ma fis = (fitssa3 + fi-1/25)/2 averaging operator 


Note that 
On = Zute de 
In a similar fashion we can define the corresponding operators in y. 
In the following table we collected some of the common approximations for 
the first derivative. 
Table I. First Derivative Approximations 


Finite Difference Order 
(See Chapter IV Section A 3) 








assis O(hz) 
Veli O(hz) 
sa behe O(n?) 
Spr (—Bhy + Afeens — Signi) = po(Ae — SAD O(h2) 
Sa Bilis — Minty + fieas) = 5-(Vet SVD O(n?) 
(hobs — SHteb)fi O(h’) 

mel oh’) 


Qheit 1§2 


The compact fourth, order three point scheme deserves some explanation. Let 


f, be v, then the method ts to be interpreted as 


a 
5p oe sii 


I 
(1 + gor) = 9 a 
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or 


l ge 
g (vit + 40,5 + Vi-1;) ae yea ‘ 


e e ° e sad s s 5) 
This is an implicit formula for the derivative st at (z;,y;). The vi; can be computed 
from the f;; by solving a tridiagonal system of algebraic equations. 
The most common second derivative approximations are 


LAS al fi — 2fiaiy + fita;) + O(he) 


HH 


feelij = (fi — 2fi-az + fi-23) + Ol he) 


1 
ta pa bedi + O(h2) 
1 Pfs | 
feelig = SS + O(hS). 
AE + 3562 


Remarks: 

1. The order of a scheme is given for a uniform mesh. 

2. Tables for difference approximations using more than three points and 
approximations of mixed derivatives are given in Anderson, Tannehill and Pletcher 


(see [Ref. 14]. 
3. Difference Representations of PDEs 


a. Truncation Error 
The difference approximations for the derivatives can be expanded in 
Taylor series. The truncation error (T.E.) is the difference between the partial deriva- 


tive and its finite difference representation. For example 


| Fins i- Sig h, 
Ie —~ —A, fi; = fr — 4 = fy a at 
Wy : 3 h, ij 2! 

We use O(h,) which means that the truncation error satisfies |JT. E.| < A h,| for 


h, — 0, sufficiently small, where K is a positive real constant. Note that O(h,) does 


not tell us the exact size of the truncation error. If another approximation has a 


23 














truncation error of O(h2), we might expect that this would be smaller only if the 
mesh is sufficiently fine. 
We define the order of a method as the lowest: power of the mesh size in 
the truncation error. Thus Table 1 gives first through fourth order approximations. 
The truncation error for a finite difference approximation of a given 


PDE is defined as the difference between the two. For example, if we approximate 


the advection equation 


by centered differences 


Fijs — Pij-1 ‘i pit — Pinas 


2At 2Az 


= 0 


then the truncation error is 

EA oe ie na eet sa ar 
ot Ox 2At 2Az 

1 OF 1 OOF 


ase 2 1 
= ~gAt ae cp Ara — higher powers of At and Az 


Todo, =( 


. We will write 
T.E. = O(At?, Az’). 
b. Consistency 
A difference equation is said to be consistent or compatible with the 
partial differential equation when it approaches the latter as the mesh sizes approaches 


zero. This is equivalent to 


T.E.—0 as meshsizes - 0. 


Cc. Stability 

A numerical scheme is called stable if errors from any source (e. g. trun- 
cation, round-off, errors in measurements) are not permitted to grow as the calculation 
proceeds. Richtmeyer and Morton give a less stringent definition of stability (see [Ref. 
15]). A scheme is stable if its solution remains a uniformly bounded function of the 


initial state for all sufficiently small At. 
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d. Convergence 

A scheme is called convergent if the solution to the finite difference 
equation approaches the exact solution to the PDE with the same initial and boundary 
conditions as the mesh sizes apporach zero. Lax has proved that under appropriate 
conditions a consistent scheme is convergent if and only if it is stable. 

The Lax Equivalence Theorem states that given a properly posed linear 
initial value problem and a finite difference approximation to it that satisfies the con- 
sistency condition, stability is the necessary and sufficient condition for convergence 


(see [Ref. 15]). 


4. Further Examples 
Given a PDE and a finite difference mesh one can use any of the following 
procedures to develop a finite difference scheme. 
a. Tables 
b. Taylor series expansions 
c. polynomial fitting 
d. integral methods 
e. control volume techniques. 
One may get the same scheme by using different approaches. As an example 


ee 0 
for procedure b we develop a three point second order approximation for a on a 
z 


Of 
nonuniform mesh. —— at point O can be written as a linear combination of values of 


Ox 
f at A,O, and B, 
0 
5 = C1 f(A) + C2of(O) + Caf(B). 
LO 





Figure 1. Nonuniform mesh 


We use Taylor series to expand f(A) and f(B) about the point O, 


Az? 


5 f'(O)+--- 


f(A) = (0 ~ Az) = (0) - Aes'(0) + = f"(0)- 
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: i ns 
f(B) = f(0 + ac) = f(0) + aAzf'(0) + = F"(0) + ——F'"(0) + - 
Thus 
Of) _ Of PO A cae Oak 
Se| = (Ci Ca + Ca)f(0) + (20s - miles | : (C1 + 0°05) S| 
L 
FCs OF Bala 
This yields the following system of equations 
Cy, +C,+ C3 =0 
Me. tulips 
1 U3 = ie 
C; - a’ C3 ==( 
The solution 1s 
a a— 1 1 
ie (a+ 1)Az’ ae aAx ’ ae a(a + 1)Az 


and thus 


Of _ —a°® f(A) +(e? —1)f(O) + f(B) 20° f 
aa aaAs + gars [to 


Note that if the grid is uniform then a = 1 and this becomes the familiar centered 


difference. 


5. Irregular Mesh 

It is more convenient to use a uniform mesh and it is more accurate in some 
cases. However. in many cases this is not possible due to boundaries which do not 
coincide with the mesh. In this case several possible cures are given in [Ref. 14]. The 
most accurate of these is a development of a finite difference approximation which is 


valid even when the mesh is nonuniform. It can be shown that 


Ae ETT iii 


Similar formula for u,,. Note that for a = 1 one obtains the centered difference 
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Figure 2. Irregular mesh near curved boundary 


approximation. 
Due to the need to refine the mesh in some of the domain to maintain the 
accuracy one is advised to use a coordinate transformation which is covered in Neta 


(see [Ref. 13]). 


6. Stability 
The problem of stability is very important in numerical analysis. There are 
two methods for checking the stability of linear difference equations. The first one 
is referred to as Fourier or von Neumann and assumes the boundary conditions are 
periodic. The second one is called the matrix method and takes care of contributions 
to the error from the boundary. | 
a. von Neumann analysis 


Suppose we solve the advection equation 


by Lax method (see [Ref. 13]) 


] At 
Fyjay = = (Fiang + Fi-nj) - ce 


9 DA gil _ Fig). 


2¢ 








If a term (a single term of Fourier and thus the linearity assumption) 
F;; = ett ptkmaz 


is substituted into the difference equation, one obtains the amplification factor 


e* = cos 3 —ivsin B 
where 
At 
Vy = c— Courant number 
Ax 
O = kyle: 
The stability requirement is 
le™*| zm ] 
and implies 
| <1. 


This is called the Courant-Friedrichs-Lewy (CFL) condition. 
b. Matriz method 
Suppose again we solve the advection equation using Lax method but now we assume 


periodic boundary conditions, 1. e. 
Faia = ie 


The system of equations obtained is 


ia — AL, 
where 
Pty 
Li 
i 
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2 : 2 
l+yp 0 l-yp 
Z 2 
A=] 0 0 
1— 
0 0 E 
2 
l-vp 0 l+vp 0 
2 2 


It is clear that the eigenvalues of A are 


Dit» ree re ; 
Aj = cos —(j —1) +wsin— (j — 1), es ee oe 


Clearly the stability of the method depends on 
(A) <1. 


Note that one obtains the same condition in this case. The two methods yield identical 
results for periodic boundary condition. It can be shown that this is not the case in 
general. See work by Hirt (see [Ref. 16]), Warming and Hyett (see [Ref. 17|) and 
Richtmeyer and Morton (see [Ref. 15]). We will explore the stability of the shallow 


water equations in more detail using Fourier techniques in the next chapter. 


7. Example: Shallow Water. Equations in 2D 

Arakawa and Lamb ([Ref. 18]) have investigated a finite difference scheme for 
the nonlinear shallow water equations using square and staggered square grids. We 
will show two of their examples, one of a square unstaggered grid, and one example 
of a staggered square grid. For more information on finite difference schemes, refer 


to their work. [Ref. 18] 
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a. Square Grid 


Recalling the linearized version of the shallow water equations 


Oh Ou Ov 
att (Se +5 | = 0, 


Arakawa and Lamb (see [Ref. 18]) use the following finite difference approximation 


for the unstaggered square grid. 


Ou g ae 

Ov IT \y — 

OL + fu+ “7 (ou) = 0 

OWE TT fees. ak 

ay tg (eu)? + (Gv)") = 0 


as ] ee 
where (6,h)* = 5 (hits — hj-1), in our notation this is p,dzh. 


b. Staggered Square Grid 

In simulation, according to Arakawa and Lamb (see [Ref. 18]), a stag- 
gered grid approach is best. The reason that a staggered square grid is better for 
modelling shallow water flow than an unstaggered square grid is that for the same 
time step, the staggered grid will give higher levels of accuracy. Now let us look at 


Arakawa and Lamb's example of a staggered square grid. Recalling the shallow water 
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equations (2.10). 


Ou -- Ou |. Ou ‘i _- — oh 

at Ax "Oy ~ ~I Oe 

Ov Ov Ov Oh 

AE re ge = 95, (4.1) 
Oh 


I 
a 


au = (uh) + = pe) 
Notice that we substituted h bai (h — a ian ignoring the bottom topography. 
For the staggered grid, following Arakawa and Lamb’s example (see [Ref. 18]), multi- 
ply the first of (4.1) by A and add it to the last of (4.1) multiplied by u, also multiply 
the second of (4.1) by h and add it to the last of (4.1) multiplied by v. This gives us 
hu, + huu, + hou, — fho + ghh, + uh; + u(hu), + u(hv)y = 0 
hy, + huv, + hvv, + fhu t+ ghhy + vhi + v(hu), + v(hv), = 0. 
This reduces to 
(uh): + (huu), + (hvu)y — fhu + ghh, = 0 
(vh); + (huv), + (huv)y + fhut+ ghh, = 0. 
This is another useful form of the momentum equations. Now multiply the first of 


2 
(4.1) by uh and add it to the last of (4.1) multiplied by = also multiply the second 
2 
of (4.1) by vh and add it to the last of (4.1) multiplied by 7 This gives us 
uhu, + hu?u, + huvu, — fhuv + ghuh, + Sh ae 5 (he (ho), = 0 
52 


vhy, + hvuv, + hv’vy + fhou+ ghvhy + +h vig (hue + 5 (hv)y = 0. 


This now reduces to 


2 2 2 
(5) + ( + (ivi — fhuv + guhh, = 0 
2/, 2). 2 , 


v? vy? v? 
("5 + (iw) + (no + fhuv + guvhh, = 0. 
2/, 2). 2} 
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This gives the equation for the time change of kinetic energy according to Arakawa 


and Lamb (see [Ref. 18]). If we multiply the last of (4.1) by gh, we get 
gh [hi + (hu)e + (hv)y] = 0 


Or 


2 


Note that the Coriolis force does not contribute to the change in total kinetic energy. 


(0% + (gh?u), + (gh?v)y — gh(uhe + vhy) = 0. (4.3) 


Also note that the sum of terms in (4.2) and (4.3) is zero. It then follows that we 
have conservation of energy which can then be used in the construction of the finite 
difference scheme. 

Arakawa and Lamb’s choice for the differencing of the continuity equa- 
tion (see [Ref. 18]) was decided in an effort to keep it as simple as possible. At h 


points, the last equation of (4.1) can be represented as 
dz Fist a F3-4j + Gi j4t = Gi;-1| = 0) 


where 
Fina = Uh ula 
Gi j4t = qh?v); 54.2 
are the mass fluxes and are defined at u and v points, respectively. This is semidiscrete 
therefore the time change terms will be left in differential form. 
For the finite difference scheme, the total kinetic energy 1s conserved 


during the inertial process. Therefore the terms 
(uh), + (huu), + (hvu), 
can be written as follows. 


a] * l nee oe m=(u)—r' >{u : 
Hu). s + aa [Se(F' at) + (GMa) + (FM a") + dy (GM uy] (4.4) 


tJ 
We will use Arakawa and Lamb’s reasoning to define H™, F™, G™, FM and GM) 


in the next few paragraphs (see [Ref. 18]). (2,7) are used as the indices and are 
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chosen in such a way that they satisfy 


O ogy. 1 = - 
= He 4 = 6. FO 4 6,G + 5 FO 4 5yGO] ~ 0) (4.5) 


Now multiply (4.5) by u;; and subtract it from (4.4), we will get 


(u) Oui, j ! Pluk 7e + CWS ay E(u a! > (u ‘ 
Hist A $,u® + GH dyuy + FM Su” + GM dius I. (4.6) 
Take (4.6) and multiply it by u;; and add it to (4.5) which is multipied by aut and a 


finite difference analog of the first three terms of the first equation of (4.2) is obtained. 
The finite difference analog is written the same as Arakawa and Lamb wrote it and 


is as follows (see [Ref. 18]). 


| 
IN 6) 92 a OO Oy a) SY 
= (H 9 U ig 92 FO juiauiag Fj i Uin1 jUing 
2 re an OO, tera tag 
a G51 Ui jUij+ — G; a Ui,j—1Ui,j 
E(u) plu 
238 2 


t I~ 9 I-95 


3(u) 3 (u) 
-, 1 Uj, 7Uj-14 ; ; } p— pues 
ee aera i,j Ui-1,j+1 Gia 51 Uid1,j 1Ui,; 


Notice that in this equation, the kinetic energy flux term appears twice, but with the 
opposite sign the second time. The definition of these terms is not dependent on the 


form of the equation as long as the total kinetic energy does not increase or decrease 


over the domain. 


The Coriolis term —fhv can be represented at u,,; by 
— f(hov");,5 


and fhu at v by 


i+dg+4 
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The pressure gradient terms ghh, and ghh, can be represented as g[h7d,h],; and 
g[h¥oyh];. ial respectively. This is a general momentum advection scheme for non- 
divergent flow to a scheme that will maintain conservation of total energy and also 
maintain the divergent flow. Horizontal errors due to discretization should be small 
on the planetary scale due to the large size of the system relative to the grid size. 
For a more detailed rendition of this example, refer to Arakawa and Lamb (see [Ref. 


18]). Neta and Lustman generalized this to nonuniform grid (see [Ref. 19]) . 


B. FINITE ELEMENT 


In the finite element methods (FEM), the domain is divided into subregions 
called finite elements, hence the name. The unknown function u is represented as the 
interpolating polynomial in each element. This representation is continuous along 


with its derivatives (to a certain order) in each elment. 


1. Basic Concepts 
One of the ways to formulate the finite element is via the so called weighted 


residuals method. In this method, the desired function u is replaced by a finite series 
; N 
u = S- 50>. 
j=l 


The set of functions ¢; are called basis functions. Clearly one can, not expect u” to 


satisfy the partial differential equation, 


ts 

= 

| 
ee 


The residual R is defined as 


R= Lu" —f. 


In order to obtain the undetermined coefficients u;, one sets the weighted residuals 


to zero, 1.e. 


| Rw; =0 i=1,2,...,N 


34 














where the weights w; may be chosen in various ways. If the weights are chosen to be 
the same as the basis functions, one obtains Galerkin’s method. For collocation, the 
weights are Dirac delta functions. In the subdomain method one uses the character- 
istic function of each subdomain as the weight, i.e. w; = 1 in the subdomain 2); and 


zero elsewhere. 


Example 


Given the equation 

2 
iid +u+z2=0 0<2x<1 
dx” 


with homogeneous boundary conditions 


we suppose that we can approximate u(x) by 


u" = ay¢1() + dodo(z) 
and let the basis functions be 


d, = x(1 — 2) 
(4.7) 
dé, = z*(1— 2). 


Notice that each basis function satisfies the boundary conditions and therefore u” 


satisfies those also. The residual is 


Rah 
a +u'+a=(-242-2*)a,4+ (2-624 2?—2)a. +2. (4.8) 


dr? 





The collocation method yields (using x; = i, f= 5 as collocation points) the 


following two equations for the unknowns a; and ap. 


29, 38, 
1671 — 6492 = 4G 
Z i ee 
4a, +3a2 = 5 
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The approximate solution is then 


42 + 40x 


A ves 
u z(1— 2x) 517 


The subdomain method yields the following system when the interval is sub- 


divided into two. 


[P Rac = 0 
0 

1 
| Rdz — 0 

1 | 
eH es, 

12° 192° 8 

11 929 3 

12° 192 ~ ~8 ace 
; _ 

; ~ 517 
. _ 8 

: ~ 517 

97 + 88x 

h _ — ee eee 
u z(1— zx) 517 


For Galerkin’s method we have the following two integrals to evaluate. 


1 
| Ré,dz = 0 
0 
(4.10) 
I 
| Rede & 0 
0 
with 
3 a + 3 i 
io 10° 12 
3 8 
10.7 105°” 20 


36 








where 


7 
“1 = 369 
ie ee ae 
7 Ad 

ao. % 
h —— — —_——. a 
w= ll —2)\se5 + a) 


The accuracy of the approximate solution depends on the choice of the basis 


functions and the method used. 


2. Weak Formulation 
The fundamental integral statement of the finite element methods can be in- 
terpreted as a combination of a weighted residual and a process of integration that 
reduces the order of continuity required. Going back to the previous example 
d*u 


aaa ee = 0 


u(0) = u(1) =i 
(4.11) 
u” = ay¢1(x) + a2¢2(z) 
1 | d?y 
| Se tut olde = P22 


Integration by parts yields 


1 


» dudd; d 
[ZZ twit vdide + Ho, 
0 


dr dz dx = eae 





It is in this weak form that one substitutes the approximate solution. Notice that 


¢;(0) = ,(1) = 0 and thus the boundary term vanishes. 


3. Choice of Basis Functions 
The accuracy of the methods of weighted residuals depends mainly on the 


choice of basis functions. In the previous example, we have chosen polynomials. We 
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now introduce several possible piecewise polynomials. The simplest polynomial is 
piecewise constants. 


a. Piecewise Constant 


a re(ti_1, 2i41) 


P(x) = 
0 elsewhere 
b. Piecewise Linear 
L— Ly-] 
————— 
Li Li-1 
Listy 
@(t)=4 “— 2<24< Liat 
Lit. — L; 
0 elsewhere 
Cc. Piecewise Quadratic 


For higher order, it is easier to introduce a nondimensional coordinate 
€ where (-1 < € < 1). This choice facilitates numerical integration by Gaussian 


quadratures. 


o_1(€) —3€(1 — £) 


bo(€) 1 ~ £’ 


o1(€) $€(1 + €) 


The subscript indicates the point € at which @ is one. © is zero at the 


l 


other two points as shown in the following figure. 
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Figure 3. Piecewise quadratic basis functions 


d. Piecewise Cubic 
It is easy to see that piecewise cubic basis functions are given in the 


nondimensional coordinate & by 


@1(€) = ye(1 — €)(9¢? ~ 1) 


AS 
An 
ti 
wT 
— 
| 


3 (3€ — 1)(€ — 1) 


(€) = —7g(3€+1)(@ - 1)» 


Oi(f) = gel + €)(96? — 1). 
The above basis functions interpolates along the element using a Lagrange third order 
polynomial. Another choice is the Hermite polynomial. In addition to continuity 
of second derivatives over the element, one has first derivative continuity between 
elements. Therefore Hermite polynomials interpolate the derivatives also. This is 
useful, for example, in flow problems where one has to differentiate the potential to 


obtain the velocity field. 


39 


Each node is identified with two basis functions. 
gr = 1(1—-€)°(€+2) — $o1(—1) = 1, dor(1) = $$1(41) = 0 
gor = (14 67(E-2) doa(1) = 1, do2(—1) = $62(£1) = 0 


1 = 11-6)°(€+1) 


do = £1486) (E-1) 








“.) 08 06 -04 =D? 0 a2 a4 a6 aa 1 


Figure 4. Piecewise cubic basis functions 


4. 2D Basis Functions 
It is very easy to extend the method of weighted residuals to higher dimensions. 
In this section we discuss several possibilities of basis functions for rectangular and 
triangular elements. We close by discussing isoparametric finite elements (irregular 
quadrilateral elements). | 
a. Rectangles 
As in one dimension, one can present Lagrangian and Hermitian basis 


functions. The Lagrangian basis functions are obtained when a product of two one 
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dimensional basis functions are used. We will use a quadratic for this example. 


iS ClO) (=) 
¢2 = 41 — @)n(1— 7) (0, -1) 
$3 = =F reo(l—H) (1,.=4) 
bs = 21-90-14) (-1,0) 
g = (1-€)(1—7’) (0, 0) 
gd = 76(1+é)(1—7’) (1,0) 
dr = —7&01-OH+n)  (-1,1) 
dg = $(1-€)n(14+n) (0,1) 


gd = 7én(1t+é)(1 +7) (1,1) 


Figure 5. Nodes associated with a rectangular element 
The 9 nodes associated with the rectangular element are shown in the 
figure above. ¢; is the basis function having a value of one at the point listed next to 


its definition and zero at the other 8 nodes. 
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It is easier to write the basis functions by distinguishing between corner 


nodes, side nodes, and interior nodes. 


Quadratic: 


Corner node 4€6(1 + €6:)nni(1 + nm) 
Side node, &=0 son(l1+nn)(1 —€) 
Side node, 7, =0 $€6 (14+ 66)(1—1”) 


Interior node (1 — €*)(1 — n”) 
Linear: 
1 
q\ + ni)(l + €&) 


Another possibility, which preserves the properties of the Lagrangian 
basis functions, but eliminates the interior nodes are called serendipity basis functions. 


For example, in the case of the quadratic: 


Corner node a(1 + €&:)(1 + n)(€& + 9: — 1) 
Side node, € =0 1(1 — €*)(1 + nm) 
Side node, hh = (0) +(1 = n7)(1 at E€;) 


The cubic and Hermite cubic can be found for example in Lapidus and 
Pinder [Ref. 20]. 

b. Triangles 

The triangular element is the most well known. It allows more accurate 
presentation of an irregular domain than the rectangular element. A natural repre- 


sentation, called area coordinates system was introduced to simplify the integration. 
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{Xz¥2) 


Figure 6. Area coordinates 


Each coordinate L; is defined in terms of the triangle sub area whose 


base is the line L; = 0 and the vertex is the point 1; = 1, thus 


where A; are the areas of the subtriangles and A is the area of the whole element. 
Clearly A = A; + Az + A3. Also clear is the L; is unity at node 7 and zero at the 
other two nodes. For example L; = 1 if P(x, y) coincides with the vertex (x,y) and 
it is zero if P coincides with either of the other two vertices. 


Note: 


i+ l,+l3=1 


m,!m!mz! 


Ly LY? Ls drdy = 2A— 
h, oe ee oe gE er ar 


where m;, m2, M3 are nonnegative integers. 


] 
a A, er ey = Za (tit? — i41)(Lj42 — Tj+41) 


] 
sp dady = — (yin — Yi 41 — Y; 
OE Oe Y aA +1 — Yit2)(Yj41 — Yt2) 
where 7 and 7 are cyclic for ] < 7,7 < 3. 
In many cases one requires or prefers numerical integration. We will 
not discuss this matter, but the reader is refered to Connor and Will (see [Ref. 21]), 


Brebbia and Connor (see [Ref. 22]) et al. 
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For quadratic basis functions, one requires three more nodes per ele- 
ment and these are taken as midpoints of each side. The six basis functions in area 


coordinates are 


403 — 1, 

AL, Ly 

213 — Le 

AL. Ls 

2132 — L3 

AL3Ly. 
Several possible higher order triangular elements can be found in Fe- 

lippa [Ref. 23], tapas and Pinder [Ref. 20], and Zienkiewicz [Ref. 24] and others. 

Cc. Isoparametric 
In the Isoparametric finite elements, all irregularly shaped elements are 


mapped onto regular elements in order to help with the integration. 


Examples 


1. Any quadrilateral with straight sides can be mapped onto a square 
whose vertices are at (+1,+1). 


2. Any quadrilateral with curved sides can be mapped onto the same 


square. 
3. Any quadratic triangle with curved sides can be mapped onto an 


equilateral triangle. 


It can be shown the the mapping from the irregular shape in z, y coor- 


dinate to &,7) is given by 
ct = ye xi; (€, 1) 


 ViGil€, 7) 


cm 
| 
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where (z;,y;) are the vertices of the quadrilateral and ¢;(&,7) are the basis functions 
of the appropriate order, (in the first case the ¢; are linear.) 


The evaluation of integrals require the Jacobian of transformation J. 


Ox Oy 

dé a€ 
J<T- 

Ox Oy 

dn On 


Isoparametric elements of higher order can be found in Lapidus and 


Pinder [Ref. 20], Zienkiewicz [Ref. 24] and others. 


5. 3D Basis Functions 


The extension to three dimensions is straight forward if one has hexahedral 
elements either Lagrangian or serendidpity (no nodes in the interior). One can also 
use tetrahedral or pentahedral elements. 

8 7 


ff 


3 





2 





73 


Figure 7. Three Dimensional Elements 


Tables for serendipity type can be found in Zienkiewicz (Ref. 24]. Verge [Ref. 


25] et al. Area coordinates are replaced by volume coordinates. 
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6. Example: Shallow Water Equations in 2D 
In this section, we give an example of solving the shallow water equation using 
finite elements. We will ignore the bottom topography. The shallow water equations 


in vorticity divergence formulation are (as one may recall from Chapter 2) 


+ (uQ), + (v@), at 
D,+ 9V?h + V7K —(vQ)c+(uQ)y = 0 
hy + [uh], + [vA], = <0; 


where fh is the geopotential height, u the east/west component of the wind, v the 
north/south component of the wind and f is still the Coriolis parameter. Hinsman 
(see [Ref. 26]) shows that the equations can be discretized using finite elements as 


follows: 
=f (ESE + SESE + State) =I [ge (wo GE) 
af -. ((uQ),¥5 - 5 ) at Pe bat (( bu) 5 Oa 2 + (ou), te 


1 OV; OV; 
anit [5 (wit a0 +(e an Fe) u— f fous Win, 








V; 








o( anit 


Ow; OV; OV; 4, OV; OV; OV; av; av, 
| “4 Ms 4 OY, BHO) __ fF (yqy,2¥) y — f (v1.2) v 


ee i ere ee eee 


Ot Ox Ox Ot Oy Oy Ot Ox Ox Ot Oy Oy 





= | Cae OV; OV; Ox; OV; OV; 06; OV; OV; — A; AV; me 


Fe eg oR N 
[& (wom -w Z| u- fF (twas — BFE) vf fousviv, 
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Divergence (V?x) = + * and = = —uf along the boundary walls. For further 
look at this example see [Ref. 26]. We have used V, for the basis functions in order 
not to confuse it with the geopotential height. Note that (4.13) and (4.14) are time 
dependent. Also note that Hinsman is using a leap-frog time discretization to obtain 
each next time step. Leap-frog requires an additional starting value obtained by the 


Matsuno scheme. (See [Ref. 5].) 
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V. STABILITY ANALYSIS 


Before getting into the stability analysis of the shallow water equations, this 
paper shall review some basics about integral transforms by focusing on the Fourier 
transforms. Integral transforms are applied in mathematical type operations by tak- 
ing an equation which is unsolvable (or really hard to solve) in its original form and 
transforming it to an equation that is algebraically solvable. Generally, the harder 
the original problem is to solve, the harder the inverse transform is to find. This 
chapter will not go into the nitty-gritty of integral transforms, but will offer some 
insight into why and how they work before applying transforms to the problem at 
hand. For more information on integral transforms see Miles [Ref. 27]. Simply put, 
an integral transform is of the form F(p) = f° K(p,x)f(ax)dz where the function to 
be transformed is f(z), and K(p,z) is the transforming function which is known as 


the kernel of the transform. 


A. ONE DIMENSION 
1. Fourier Transform 


The Fourier Transform in one dimension is 


i(k) = F(f(a)) = fe f(a)de. 


— 


Also, recall that the inverse Fourier Transform is 


f(z) = FO (f(k)) = (—) f. j(kei*ak. 


2. Shallow Water Equations 

Schoenstadt (see [Ref. 28]) studied the effect of replacing the spatial deriva- 
tives in a dispersive wave equation by using Fourier transform techniques. We will 
repeat some of his work here, but for a more in depth understanding, the reader 


should consult the original text. Let us consider the effect of semi-discretization on 
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the solution of a shallow water type model. The shallow water equations in a one- 


dimensional linearized form with no mean flow, as one may recall from Chapter 3, 


are 
u:— fut+gh; = 0 
v+ fu = 0 (5.1) 


Now if we use the Fourier transform on these equations, (hats (*) symbolize Fourier 
transforms), we arrive at 

i, = fb—ikgh 

‘ _ fii 

hy = —tkHt 


= 
| 


The initial conditions are 


0 
dp = 6(k, 0) = f[R v(z, Oe“ dx 
ho = h(k, 0) = f%, h(x, 0)e~** dx 
The transformed shallow water equations are a coupled set of constant coefficient 


ordinary differential equations which can be solved fairly easily (see [Ref. 28]) and 


are 
7 M9 .- Oo : 
“= a2 pit = —3 pw ivt 
Vy Vy 
i . wpa wW t Q3 iv 
6 = tkga,+ be -_ Jind (5.2) 
Vy y 
, kHa2 ; kHaz _. 
h pe fa _ > 2 eit + > 3 Tit 
Vy Vv 
where 
v? = f?4k?gH = f2(1 + dk?) 
» = vg 
f 
; . , ; bts ns ; ee ae 
The a,’s are picked to satisfy the initial conditions. Recalling that sins = aaa 
. ; 1 
Cape.” , ; ins sc 
and cos r = ——5—— we can rewrite (5.2) by using the intial conditions to solve for 


a_a 
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a; and then simplifying to get 


i(k, t) 


o(k, t) 


h(k, t) 





The steady state solutions are 


or to rewrite 


By noting that 


tip cosvt + Jo sin vt — akghto sin vt 
y V 
P 2 2 
_ ies sinvt + aes F cos u| Vo + saat {1 — cos vt} ho 
kt ikH ROL , 
eee ut — 2 1 — cost} i + a as ecavi ho. 
V y? 
u,(k) = 0 
; k’gH, — ikgf ; 
d5(k) = 2 Vo ey ho 
s kA f . ie 
h,(k) —" = 2 Vo — “a0 
u;(k) = 0 
‘ . , f? frkg, 
d.(k) = tot+ =) Fe — Vo 
, ~ AA f*, fikg,s , 
s(k) = hot Paro (“Fig = i] 
=i] _ spp 2X 2A f? 
ee te are y2 




















ol 





and using the convolution theorem, Schoenstadt (see [Ref. 28]) inverted the steady 


state transformed equations to yield 


u;(z) = 0 


v(x) = v(z,0)+ = [ a as (f(s ,0) — o(s,0) ds 





He ah 
h,(x) = h(x, 0) — sail. sen(z at (25 oe 


Note that in the one dimensional case, u(k,0) does not contribute to the steady state 


(s,0) — v(s, 0) ds. 


solutions. Schoenstadt (see [Ref. 28]) goes on to explore the transient parts of the 
model, and these can be studied in more detail in his 1977 paper. 

Schoenstadt (see [Ref. 29]) also analyzed a difference scheme for one dimension 
shallow water equations that is similar to the one we will present later in this chapter 


for two dimensions. 


B. TWO DIMENSIONS 
1. Fourier Transform 


Recall that the Fourier Transform in 2 dimensions is 


f(wi,w2) = F(f(e,y)) = ff f(e,yesonM dedy, 


and the inverse Fourier Transform 1s 


f(x,y) = (= -) f mi (w 1, we)e ial ad ae 


Also recal] several elementary properties such as: 


Fit) = SFU 
Fife) = iwF(f) 
Fy) = iuF(f) 
FIV) = iF({) 
F(V*f) = —w FUN) 
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We will use several of these properties later in this chapter when dealing with the 
continuous case and also the semi-discrete case. The convolution theorem will also 


be used in this section and, in two dimensions, it is 


f(x,y) = Fo" (g(k, 1) ; h(k.1)) = [fF  9(0,vo)h(a — Lo, Y — yo)dxodyo. 


2. Continuous Case 


Recall the linearized shallow water equations from Chapter 3. 


Ou Oh 

a T° + 99, =" 

0 Oh 

Ge ag, (5.3) 
Ad 2 bs Ga er 

Ot Ox Oy} 


Following the work of Neta (see [Ref. 2]), we will go through the reasoning 
done by Neta and by Neta and DeVito (see (Ref. 7]) to arrive at their solution. Let 


us take the first equation of (5.1) and transform it using Fourier, 


0 h . 
* (Gi fosge) = [I (Gi sot age) emer 


= = / ~ | 7 ue +) dady — f | ° / ~ ve Ett!¥) dady + g | 7 ‘ . oe ete) dy 
—00 J —00 —co J —00 —ood—co OL 


= fit ighh = 0 


Similarly 
Ov Oh Ov er, 
F (H+ fut oo" | = 3, t fut igh = 0. 


Oh du Av Oho. 
F (Spt n (Se 4S)) =F sinha + imo =o. 


ay) 








This is the reasoning that Neta and De Vito (see [Ref. 7]) used to show that the 


Fourier transform of (5.1) is 


U us uo Uo 

1 ] , 
6 {=| 8, 1+ vat By | cosvt + 3B vo | sinvt (5.4) 
h h, ho ho 


where A and B are given by 


ft +ghk? gHkli wg fl 
A= | gHkl f?+ghl? —igfk 
—iHfl tH fk gH(k? +1?) 
and 
0 f —igk 
B=| —f 0 —igl 
~tHk ~—tHl 0 
Note that if you reduce the terms to one dimension, the results are the same as shown 


by Schoenstadt in the one dimensional case. The steady state solution for (5.4) is 


Us Uo 
1 
vs mar Vo 
h, ho 
where | 
ghl?  =—gHkl —igfl 
C=] —gHkl ghk? — igfk (5.5) 


tHfl  —-iHfk f? 
up, Uo and ho are the Fourier transforms of the initial conditions. The frequency v 
is given by v = f,/1 + A2(k? + [*) where A is the Rossby radius of deformation (i.e. 
hte 
Following Neta (see [Ref. 30]), we shall now take the inverse Fourier transform 


of (5.5) to obtain the steady state solution. We will follow the arrangement of the 


o4 

















variables used there. They are 
D = iktip + ilig 


= tltig — tk 


I> 


f 


d, —_ i2 kh ae Vo 
d, = Uh + tio 


where the steady state solution is 


tile = ty + S(NAkD — d,) 
2 
v, =Uvot L ip + dz) 
VY 
i a Hf? jie e , d. 


We can now invert these equations using the convolution theorem and the following 








integrals. 
i 2x Ct 
4! epee Oddb = Io(rp) 
1 re Jo(rp) 1 r 
oe | T+ Az 2? = oe a (5.6) 
2 
= ie 1 Jax? + y? 
- (= = ako d 


Jy and Ao are Bessel functions of order zero. Using these integrals, the inverse 


obtained from them, and the convolution theorem, we have the steady state solutions. 


a3) 











1 re re € 
Us(2, y) uo(z, y) =~ | fs Kolae, T))D(o, T, 0)dodr 
1 CO CO 7 
-sa/_ | Kolao,7))a(o, T, 0)dodr, 


1 po poo 8 
v(2y) = volzy)—s-f f <-Kola(o,7))D(o,7,0)dodr 


] CO CO 
to—5] f Kola(o,7))delo, T, 0)dedr, 


and 


h,(x,y) 


ho(z,y) — saul. [ (3 Ko(q(c,r))d,(o, 7,0) 


ee ae (0,7) )dz(o, 7 0) dodr 


Oz 
where 
(x—o)° +(y—T) 
q(o, T) _ 
O O 
D(c,7,0) = — +3 
dh 
d,(o,7,0) = — ~ u% 
ol 
d,(o,7,0) = en Up. 


3. Semi-discrete Equations 
Here we will present the semi-discrete shallow water system and then show the 
inverse Fourier transform for certain selections of a, 3, and y’s. Semi-discrete means 


that time is not discretized. In the semi-discrete case, the shallow water equations 
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are 


Ou St ae ee 

az —Pfotignh = 0 

ax + Bfitignh = 0 (5.7) 
Oh oo 

an +72H(yt+72b) = 0 


where a, J, 71, and 72 depend on the scheme used (see tables II and III) from Neta 
[Ref. 30]. 


Table II. Filter weights for second and fourth order finite differences 

















Scheme a B v1 % 
sin kd sin ld 
A2 ] 1] 
d d 
kd... Id E kd 
sin = cos = sin = cos = 
B2 ] ] a 2 x 
2 
. kd ' 
C2 COs — Cos i ee sin > 
: 2 d 
2 2 
D2 ] kd ld sin kd cos = ig sin ld cos = i 
fd cos — nee 
ST cos > ar ca ; 
A4 l 1 8sin kd — sin 2kd 8sin ld — sin 2ld 
~ ca * ees 
BA I (csin 3 + 27sin f)cosF (sin St + 27 sin Z) cos 9 
12d 12d 
C4 cos — cos i — sin ois + 27 sin — sin a + 27 sin = ‘ 
2 2 12d 12d 
DA hed ld (8 sin kd — sin 2hd) cos “g (s sin add — sin 2ld) cos xg 
; cos ~~ Cos — SSA near Oe Me ee ae re eh ee oor 
= 6d 6d 


of 





Table III. Filter weights for finite elements 


Scheme a= B “1 . 
pet Bt coskd + 2cos 7 cos ld i cos Ba ld 
6 3 


(2 + cos kd)(2 + cos Id) (2 + cos ld)(sin kd) (2 + cos kd)(sin ld) 


= 9 3d 3d 


The two schemes for finite elements in Table III are for triangles (Finite 
Element-Triangle or FET) and rectangles (FER). 

The Fourier transform of the steady state solution is similar to the steady state 
solution of the continuous case and is written as follows with obvious changes from 


the continuous case. 


Us uo 
1 A 
Vs Ate D Vo 
h, ho 
where 
ghy3 —gHy172 —igf By 
Cpo=| -—gHy172 ghy? igf By: (5.8) 
itHfBy. —-iHfBy, Bf? 
and 
avy = Bf V/1+ADB(¥7 +79) 
a = 


Bf 
There is a great deal of similarity between the continuous case and the discrete case 


up to this point. We can carry it further and define as before 


D(k,1,0) = inte + i726 
((k,1,0) = iyato — indo 
d,(k, l, 0) = i rho — Boo 


d,(k, l, 0) — iF ake + Buo. 
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Therefore the steady state Fourier transforms of the semi-discrete shallow water equa- 








tions are p 
fi, = to + GV inD — Bay) 
2 me : 
Vs = Uo + aie (A wy2D ea Bd, ) (5.9) 


h, = ho + 5 lind. + i72dy). 
Notice that the Rossby radius of deformation (A) was used instead of the its discrete 
analog (Ap). 
| a. Finding the Steady State Solutions for Scheme A2 


Let 


et(krt+ly) 
Im(z,y) = F+(55)- [- of ee ED 4 MED LaAE 7 hai 


For scheme A2, F—'(i414(k, L)) is 








_1/sinkd ., © £4 Oq(xo, yo) 1 
se or ra) ie 3 Bx 2d0\Y ~ Yo)dtodyo 


2dJ-d Ox 
ro =jla(d,y) ~ o(-4,y) 


and Fo! (ty29(k, 1)) is 





_,,sinld Oq(zo, 
FS igk)) = [- [Mate 1 se — xo)dyodzc 


ld 
_ ae (x Yo) 
7 al’, 
= = ,d) — d 
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It follows in the semi-discrete case of the shallow water equations, by using the con- 


volution theorem the steady state solution for scheme A2 is 


u,(r,y) = uo(z,y) + so oe htt [uo(d, 7,0) — uo(—d, 7, 0) 
1\2 
—I,,(z — o0,y—T) (=) [{vo(d, d) — vo(d, —d)} — {vo(—d, d) — vo(—d, -<))| dodt 
-[ f m(x — 0,y —T)dy(o, 7, 0)dodr 
V(r, y) = vo(z,y) + $f f. pei ee a os. 57 lvole, d,0) — vo(a, —d, 0)] 
Si iGgee ven (=) FONG eT Geen ee CT re ee -<)y)} dodr 
ff in(e m(t — o,y — T)d,(o,7, 0)dodr 
hay) = hola of” [” | Sale 2d Shad 1,0) — ho(—d, 7,0) 


—I,,(r — o,y — 7)[vo(d, 7,0) — vo(—d, 7, 0)] 


—1,,(2 —oO,y—- T)[uo(o, d, 0) ~~ uo(d, ~—d, 0)] 


Ol,,(r -—o.y — 
: (x y—T) 


7lholo: d,0) — ho(o, —d, 0) dodr 


Oy 
(5.10) 
where ec 
d,(z, y, 0) = f 2d [ho(d, y) =, ho(—d, y)| — Uo 
i 
d,(z,y,0) = EP [ho(z, d) — ho(x, —d)] + uo. 
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Note that d,(z,y,0) is a centered difference approximation which approaches u, as 
d approaches 0. One should be able to obtain similar results for the other schemes 
listed in the tables. 

b. Solving for In 

The reason for this thesis is to show which semi-discrete method has 
an inverse Fourier transform which can be back transformed in a closed form. By the 
Lax equivalence theorem, if we can show that J,, is bounded in some sense, then the 
solution is stable also. Numerical analysis may be used to solve this problem, but a 


simple closed formula is considered a more ideal solution. Let us look again at 


et(ketly) 
Im(2,y) = FO ll =p ee 
owd= (55) = [mene 6 
Note that the limit of [,, for any of the schemes when d -+ 0 is 
et(ke+ly) 
dkdl. 
~=| oie 
By using the identity 4 = f¢° e~**dA, we have the following triple integral, 
‘ea = ia ia eilketly) {ate (M4) ap didu, 
0 —oo/ —00 
which can be rewritten as 


L,, = | ae (| ake cilve oF a) dit. 
0 —0Oo —~0O 


Evaluate the three integrals separately in this form and it easily can be shown that 





] 
] 2 EEE ‘é 
lim | 5 Ko 


where Ao is the Bessel function of order zero from (5.6). This shows that the steady 
state solution of the semi-discrete system (5.10) approaches the corresponding solution 
(5.5) of the continuous case. 

Refering to Neta and Devito (see [Ref. 7]) we have several second and 


fourth order finite difference equations (Table II). In scheme A2 we have 
] ff t(kr+ly) 
A2 = ry a ON 5.12 
° -coJ—oo] + )2 (ee + a) ( ) 
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By noting in (5.12) that we have, for example, cosk in the numerator and that the 
denominator is positive and finite, we see that [42 does not converge. What one can 
do is look at [42 to determine if there is any useful information which can be drawn 
from it. This information then can be used to predict some aspect of the semi-discrete 
shallow water equations. Latta ([Ref. 31]) suggested the following reasoning. One 


can rewrite [,. as 
rey ey dkdl 
ln= ol [3 
aC 00 4 (sin? (kd) + sin (Id)) 


where C' = a Notice that C is a constant. Now if one lets kd + k and Id > 1, Iy2 


becomes 


* 5.13 


Letk 9 k+a7 andl —-l+4+7. This ieee e _ is infinite where I4. = e*¢"I 42 and 
Tao = e'@7 I 49. Therefore [42 is infinite when x = 2md and y = 2nd where m and n 


are integers and J42 = 0 otherwise. The solution to [42 will look like 


oo 


7. = _ _ - I 
Cd*Iaq = Se S35 Omnd(k 2m)d(1 an) = F (re aa) 


m=—oOoOn=—CoO 


Now take the inverse Fourier Transform of the middle term to get 


it 7 <= = t2mk_ 42nl I 
fs d dy Amner™e ~ Ty sin?(k) + sin*(1) 


By noting that we are only interested in the real part of [42, we can write the equation 


as 
= = 47? 


Amr, cos (2mk) cos (2nl) = ——_----_——-~. 
oon Om O08 (2mm) 08 (Pal) = Fae) + aim 
Using the Fourier Cosine series, we can write the coefficients Qmn as 


2m pln 
Onn = | Ps cos (2mk) cos (2nl) kdl. 
& +sin’(k) + sin’(l) 


Now using the trigonometric substitution sin? 6 = ae we get 
an 2m cos (2mk) cos (2nl) 
maf i, Ty _ cost) + costa) 
ot 2 





which can be written as 


m= fe ig cos (2mk) cos (2nl) ddl. 
1+C 1 _ £08 (2k) + cos (21) 
ie oi +C 
C 





Using the geometric series expansion (5 =1l+z2r+2*4+2°+...), we can write 
—2 


Am 2) as 


an [Pont cos (2mk) cos (2m) ere i (cos (2k) + cos (21)) dkdl (5.14) 


Note that cos (2mk) cos (2nl) < 1 and (cos (2k) + cos (21))” < 2”. Therefore 


Amn < 4m?) Nee <a) . 


This means that an, is bounded and it follows that although J, infinite, it is stable 


_at each lattice point. 


r=0 


a V 
Since (a+ 6)” = > ( a’~"b", (5.4) becomes 
, | 


oO C Voev 


V 27 27 
Omn = » (<5) = ; / cos (2km)cos (2k) dk | cos (2In)cos”” (21)dl. 


Using Gradshteyn and Ryshik (see [Ref. 32], p. 374), one gets the following definite 


integral. 


Io. cos (2mk)cos"(2k)dk = 


r! 


* (2m —r)(2m —r + 2)...(2m +r) som; 
_ 
(pAb (1 em") =i [2m <r and r — 2m = 21]; 
r! 


[2m <r and r—2m = 21 + 1] 
(5.15) 


(21+ 1)N(4m +248 
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where 


0 [2m —r = 2:2]; 
$=, 1 [m-—r=4i+4+ 1); 
—-1 22m-r=41-]] 


Since Qmn is the same across any 27 interval, we can write it as 


v=0 7=0 Tr =, 


) 3 (Ss) o ( g | [. cos (2km)cos" (2k)dk cos (2/n)cos’~” (21) dl. 


Noting that cosine is even, this can be rewritten as 


r=0 


CO C VY Vv V T Rr 
Omn = a, (sa) > ( ; | [ cos (2km)cos" (2k) dk f cos (2In)cos”~" (21)dl. 
By using (5.15), we note that if v is odd, then an» = 0. Also if v is even and r is odd 
then Qmn = 0. Therefore vy = (2m + 2n + 27 + 27) where r = 2m + 22. Now we note 
ifr <2mory—r < 2n, then amn = 0. From these observations of (5.15), it follows 


vis ce Poe V T r T y—r 


that 


v=0 
This can be reduced to 


‘eS ie y! 


. 7) (r — aly —r —7)'y! 


v=2m+2n r=2m 


Where the prime on the summation means to count only even indices. Therefore 


x? a0 ox, fore) Pop ft C Vv y! , 
tan = OH s- De, > \, (a) Che 7s ee 


M=—CONn=—COv=2m+2n r=2m 
or in terms of A 


Oo ore) top i 2 | 
( : = 6(k—2m)(1—2n). 


m= SS YY (gecm) Goa 


mM=—-COON=—COV=2mM4t2n r=2m 
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VI. CONCLUSIONS 


The basis of this thesis was to help the reader form a better understanding 
of the shallow water equations and also to see if an inverse Fourier transform for 
the semi-discretized system can give any information on how the system behaves as 
a function of time. The first four chapters consisted of background material that ] 
compiled as a result of studying the shallow water equations, and also in an attempt to 
give the reader a logical order to better facilitate understanding of the later material. 
As can be seen, various forms of the shallow water equations were discussed and 
derived. Most of that work was taken directly from the references and were so noted. 
I feel this thesis will stand as a means for future students to more easily come to an 
understanding of shallow water equations and help them to pursue the next step in 
the process by finding a way of looking at other finite difference and finite element 
schemes. 


I hope you enjoyed reading this thesis. 
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